#    This file is part of EAP.
#
#    EAP is free software: you can redistribute it and/or modify
#    it under the terms of the GNU Lesser General Public License as
#    published by the Free Software Foundation, either version 3 of
#    the License, or (at your option) any later version.
#
#    EAP is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#    GNU Lesser General Public License for more details.
#
#    You should have received a copy of the GNU Lesser General Public
#    License along with EAP. If not, see <http://www.gnu.org/licenses/>.

import random
import sys
import math
import copy
import matplotlib.pyplot as pyplot

from eap import base
from eap import creator
from eap import toolbox


class Landscape(object):
    """A 2-dimensional bounded graph of hotspots generated from one of Sherri's dat files"""
    def __init__(self,filename):
        lines = open(filename).readlines()
        self.xMin = float(lines[1])
        self.xMax = float(lines[2])
        self.yMin = float(lines[3])
        self.yMax = float(lines[4])
        self.hotspots = []
        for line in lines[5:]:
            coords = line.split()
            self.hotspots.append([float(i) for i in coords])

def randPoint(landscape):
    """Generate a coordinate pair located somewhere within a landscape."""
    return [random.uniform(landscape.xMin, landscape.xMax), random.uniform(landscape.yMin, landscape.yMax), 0.0]

def distance(a,b):
    """Euclidean distance from a to b."""
    xterms = (b[0] - a[0])**2
    yterms = (b[1] - a[1])**2
    return (xterms + yterms)**0.5
    
def minmax(lower, value, upper):
    """Forces a value to fit with the range lower to upper"""
    value = min(value, upper)
    value = max(value, lower)
    return value

def evalTowerPlacement(individual, landscape=None):
    """The basic fitness evaluation used by the hypothetical company."""
    fitness = 0
    for shots,dartboard in [[individual,landscape.hotspots],[landscape.hotspots,individual]]:
        for dart in shots:
            dist=landscape.xMax + landscape.yMax
            for point in dartboard: 
                if distance(dart, point) < dist:
                    dist = distance(dart, point)
            if dist > 1: fitness -= dist
            else: fitness += (1 - dist)

    return fitness,

def evalTowerPointwise(individual, landscape=None):
    """Evaluate the fitness of each point and store it."""
    fitness = 0
    remaining = set(tuple(i) for i in landscape.hotspots)
    for tower in individual:
        score = 0
        covered = set()
        dist = landscape.xMax + landscape.yMax
        for spot in remaining:
            dp = distance(tower, spot)
            if dp < 1:
                score += 1
                covered.add(spot)
            if dp < dist: dist = dp
        if score: tower[2] = score
        else: tower[2] = -math.floor(dist)
        if dist > 1: fitness -= dist
        else: fitness += score
        remaining -= covered
    return fitness,

def mutSinglePoint(individual, indpb, landscape):
    """Modify the components of a single randomly chosen point with probability indpb."""
    point = random.sample(individual,1)[0]
    
    if random.random() < indpb:
        point[0] += random.gauss(0, landscape.xMax)
        point[0] = minmax(landscape.xMin, point, landscape.xMax)

    if random.random() < indpb:
        point[1] += random.gauss(0,landscape.yMax)
        point[1] = minmax(landscape.yMin, point, landscape.yMax)

    return individual

def mutSeveralPoints(individual, indpb, landscape):
    """Choose k points, where k is from a distribution centered around the number of towers/10,
    and mutate their components with probability indpb."""
    k = int(math.floor( minmax( 0, random.gauss(len(individual)/10, len(individual)/3), len(individual))))
    points = random.sample(individual,k)
    for point in points:
        score = point[2]
        sigmaX = landscape.xMax/3
        sigmaY = landscape.yMax/3
        if score > 0:
            sigmaX /= (score**2 * sigmaX**3)
            sigmaY /= (score**2 * sigmaY**3)
        else:
            sigmaX -= score
            sigmaY -= score
        if random.random() < indpb:
            mod = random.gauss(0, sigmaX)
            #print "Modifying X=%s by %s, score is %s" % (point[0], mod, score)
            point[0] += mod
            #point[0] = minmax(landscape.xMin, point, landscape.xMax/3)
            point[0] = randPoint(landscape)[0]

        if random.random() < indpb:
            mod = random.gauss(0, sigmaY)
            #print "Modifying Y=%s by %s, score is %s" % (point[1], mod, score)
            point[1] += mod
            #point[1] = minmax(landscape.yMin, point, landscape.yMax/3)
            point[1] = randPoint(landscape)[1]

    return individual

def mutPoints(individual, indpb, landscape):
    """Mutate all points wholly with probability indpb."""
    xFactor=landscape.xMax/3
    yFactor=landscape.yMax/3
    for point in individual:
        if random.random() < indpb:
            point[0] += random.gauss(0,xFactor)
            point[1] += random.gauss(0,yFactor)
            point[0] = minmax(landscape.xMin, point, landscape.xMax)
            point[1] = minmax(landscape.yMin, point, landscape.yMax)

    return individual

def main():
    #random.seed(64)
    
    pop = tools.population()
    CXPB, MUTPB, NGEN = 0.5, 0.2, int(sys.argv[2]) 
    
    print "Start of evolution"
    
    # Evaluate the entire population
    fitnesses = map(tools.evaluate, pop)
    for ind, fit in zip(pop, fitnesses):
        ind.fitness.values = fit
    
    print "  Evaluated %i individuals" % len(pop)
    evals = len(pop)
    
    # Begin the evolution
    for g in range(NGEN):
        print "-- Generation %i --" % g
        
        # Select the next generation individuals
        offsprings = tools.select(pop, n=len(pop))
        # Clone the selected individuals
        offsprings = map(tools.clone, offsprings)
    
        # Apply crossover and mutation on the offsprings
        for child1, child2 in zip(offsprings[::2], offsprings[1::2]):
            if random.random() < CXPB:
                tools.mate(child1, child2)
                del child1.fitness.values
                del child2.fitness.values

        for mutant in offsprings:
            if random.random() < MUTPB:
                tools.mutate(mutant)
                del mutant.fitness.values
    
        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ind for ind in offsprings if not ind.fitness.valid]
        fitnesses = map(tools.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit
        
	evals += len(invalid_ind)
        print "  Evaluated %i individuals, %i evaluations so far" % (len(invalid_ind), evals)
        
        # The population is entirely replaced by the offsprings
        pop[:] = offsprings
        
        # Gather all the fitnesses in one list and print the stats
        fits = [ind.fitness.values[0] for ind in pop]
        
        length = len(pop)
        mean = sum(fits) / length
        sum2 = sum(x*x for x in fits)
        std_dev = abs(sum2 / length - mean**2)**0.5
        
        print "  Min %s" % min(fits)
        print "  Max %s" % max(fits)
        print "  Avg %s" % mean
        print "  Std %s" % std_dev
        
        if g % 25 == 0:
            outfile = open('out','a')
            outfile.write("%s\n" % evals)
            outfile.close()

        # Increase the mutation rate based on the standard deviation.
        # Problem: the "good" range of the standard deviation tends to change a
        # lot for different methods.
        # Also that I don't really know any stats.
        #if std_dev < 1.5: 
        #    tools.unregister("mutate")
        #    mutrate += 0.005
        #    tools.register("mutate", mutator, indpb=mutrate, landscape=landscape)
        #elif std_dev > 4:
        #    tools.unregister("mutate")
        #    mutrate -= 0.005
        #    tools.register("mutate", mutator, indpb=mutrate, landscape=landscape)
        #mutrate = max(0.05, min(mutrate, 0.8))
        #print "Mutation rate: %f" % mutrate
    
    print "-- End of (successful) evolution --"
    
    best_ind = toolbox.selBest(pop, 1)[0]
    #print "Best individual is %s, %s" % (best_ind, best_ind.fitness.values)
    for tower in best_ind:
        print tower
    print evalTowerPlacement(best_ind,landscape)
    print evals
    towerX = [i[0] for i in best_ind]
    towerY = [i[1] for i in best_ind]
    hotX = [i[0] for i in landscape.hotspots]
    hotY = [i[1] for i in landscape.hotspots]
    
    #pyplot.scatter(towerX,towerY,c='r')
    #pyplot.scatter(hotX,hotY)
    #pyplot.show()



if __name__ == "__main__":
    creator.create("FitnessMax", base.Fitness, weights=(1.0,))
    creator.create("Individual", list, fitness=creator.FitnessMax)
    creator.create("Population", list)
    
    landscape = Landscape(sys.argv[1])
    
    tools = toolbox.Toolbox()
    # Attribute generator
    tools.register("attr_point", randPoint, landscape)
    # Structure initializers
    tools.register("individual", creator.Individual,
        content_init=tools.attr_point, size_init=30)
    tools.register("population", creator.Population,
        content_init=tools.individual, size_init=300)
    
    # Operator registering
    tools.register("evaluate", evalTowerPointwise, landscape=landscape)
    #tools.register("evaluate", evalTowerPlacement, landscape=landscape)
    tools.register("mate", toolbox.cxTwoPoints)
    #tools.register("mutate", mutPoints, indpb=0.05, landscape=landscape)
    mutator = mutSeveralPoints
    tools.register("mutate", mutator, indpb=0.33, landscape=landscape)
    tools.register("select", toolbox.selTournament, tournsize=3)
    #tools.register("select", toolbox.selBest, n=15)
    main()
